% This file analyzes the experiments intended to show the phase transition
% [r,0,-r]

%% read all relevant data
load('./data/FigS15to20.mat');
name_PT={'025','045','06','0825','1','sqrt2'};
r=2e6;
flag = {[1,0,0],[0,1,0],[0,0,1],[1,1,0],[1,-1,0],[1,0,1],[1,0,-1],[0,1,1],[0,1,-1]};
beta=0;phi=0;

for hlp_PT=1:6    
    %% plot
    alpha_fit=linspace(0,pi/2,41);
    [rabi_estimate_save,g_munu_theory,H_theory] = QGT_rabi_freq_detune_theory(alpha_fit,0,0,rd(hlp_PT),2e6,flag);
    rabi_estimate_save=rabi_estimate_save/1e6;
    % Rabi oscillations
    figure
    cc=[0.9647    0.3294    0.0510
        0.9020         0    0.0627
        0.0941    0.5294    0.9098
        0.8941    0.2314    0.5333
        0.1020    0.5686    0.4745
        0.9882    0.7098    0.0588
        0.5686    0.7490    0.1020
        0.0549    0.2706    0.6275
        0.0980    0.5412    0.1843];
    
    % SQ
    hold on
    for hlp=1:length(flag)
        switch hlp
            case {4,6}
                plot(alpha_fit/pi*180,reshape(rabi_estimate_save(1,hlp,:),1,[]),'--','LineWidth',1.5,'Color',cc(hlp,:),'HandleVisibility','off');
            case {5,7}
                plot(alpha_fit/pi*180,reshape(rabi_estimate_save(1,hlp,:),1,[]),':','LineWidth',1.5,'Color',cc(hlp,:),'HandleVisibility','off');
            otherwise
                plot(alpha_fit/pi*180,reshape(rabi_estimate_save(1,hlp,:),1,[]),'LineWidth',1.5,'Color',cc(hlp,:),'HandleVisibility','off');
        end
    end
    
    for hlp=1:length(flag)
        switch hlp
            case {1,2,3} % diagonal terms, solid circles
                errorbar(alpha_sweep{hlp_PT}/pi*180,reshape(rabi_fit_save_PT{hlp_PT}(1,hlp,:),1,[])./reshape(E_pert_save_PT{hlp_PT}(1,hlp,:),1,[]),reshape(drabi_fit_save_PT{hlp_PT}(1,hlp,:),1,[])./reshape(E_pert_save_PT{hlp_PT}(1,hlp,:),1,[]),'o','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:),'MarkerFaceColor',cc(hlp,:));
            case {5}
                errorbar(alpha_sweep{hlp_PT}/pi*180,reshape(rabi_fit_save_PT{hlp_PT}(1,hlp,:),1,[])./reshape(E_pert_save_PT{hlp_PT}(1,hlp,:),1,[]),reshape(drabi_fit_save_PT{hlp_PT}(1,hlp,:),1,[])./reshape(E_pert_save_PT{hlp_PT}(1,hlp,:),1,[]),'s','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:),'MarkerFaceColor',cc(hlp,:));
            case {4,6}
                errorbar(alpha_sweep{hlp_PT}/pi*180,reshape(rabi_fit_save_PT{hlp_PT}(1,hlp,:),1,[])./reshape(E_pert_save_PT{hlp_PT}(1,hlp,:),1,[]),reshape(drabi_fit_save_PT{hlp_PT}(1,hlp,:),1,[])./reshape(E_pert_save_PT{hlp_PT}(1,hlp,:),1,[]),'s','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:));
            case {7,8}
                errorbar(alpha_sweep{hlp_PT}/pi*180,reshape(rabi_fit_save_PT{hlp_PT}(1,hlp,:),1,[])./reshape(E_pert_save_PT{hlp_PT}(1,hlp,:),1,[]),reshape(drabi_fit_save_PT{hlp_PT}(1,hlp,:),1,[])./reshape(E_pert_save_PT{hlp_PT}(1,hlp,:),1,[]),'^','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:));
            case {9}
                errorbar(alpha_sweep{hlp_PT}/pi*180,reshape(rabi_fit_save_PT{hlp_PT}(1,hlp,:),1,[])./reshape(E_pert_save_PT{hlp_PT}(1,hlp,:),1,[]),reshape(drabi_fit_save_PT{hlp_PT}(1,hlp,:),1,[])./reshape(E_pert_save_PT{hlp_PT}(1,hlp,:),1,[]),'^','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:),'MarkerFaceColor',cc(hlp,:));
            otherwise
                errorbar(alpha_sweep{hlp_PT}/pi*180,reshape(rabi_fit_save_PT{hlp_PT}(1,hlp,:),1,[])./reshape(E_pert_save_PT{hlp_PT}(1,hlp,:),1,[]),reshape(drabi_fit_save_PT{hlp_PT}(1,hlp,:),1,[])./reshape(E_pert_save_PT{hlp_PT}(1,hlp,:),1,[]),'o','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:));
        end
    end
    
    lgd=legend('$$\Gamma_{\alpha}$$','$$\Gamma_{\beta}$$','$$\Gamma_{\phi}$$','$$\Gamma_{\alpha\beta}$$','$$\Gamma_{\alpha\overline{\beta}}$$',...
        '$$\Gamma_{\alpha\phi}$$','$$\Gamma_{\alpha\overline{\phi}}$$','$$\Gamma_{\beta\phi}$$','$$\Gamma_{\beta\overline{\phi}}$$');
    set(lgd,'Interpreter','latex');
    set(lgd,'Location','eastoutside')
    lgd.FontSize=18;
    
    ymax=(max(max(rabi_estimate_save(1,:,:))))+0.25;
    ylim([0,ymax])
    xlim([0,90])
    xticks(0:15:90)
    xlabel('\alpha(\circ)','FontSize',18)
    ylabel('Matrix Element (MHz)','FontSize',18)
    title('SQ','FontSize',20)
    hold off
    set(gca,'FontSize',18)
    
%     saveas(gcf,['SQ_',name_PT{hlp_PT}],'epsc')
    % DQ
    figure
    hold on
    for hlp=1:length(flag)
        switch hlp
            case {4,6}
                plot(alpha_fit/pi*180,reshape(rabi_estimate_save(2,hlp,:),1,[]),'--','LineWidth',1.5,'Color',cc(hlp,:),'HandleVisibility','off');
            case {5,7}
                plot(alpha_fit/pi*180,reshape(rabi_estimate_save(2,hlp,:),1,[]),':','LineWidth',1.5,'Color',cc(hlp,:),'HandleVisibility','off');
            otherwise
                plot(alpha_fit/pi*180,reshape(rabi_estimate_save(2,hlp,:),1,[]),'LineWidth',1.5,'Color',cc(hlp,:),'HandleVisibility','off');
        end
    end
    
    for hlp=1:length(flag)
        switch hlp
            case {1,2,3} % diagonal terms, solid circles
                errorbar(alpha_sweep{hlp_PT}/pi*180,reshape(rabi_fit_save_PT{hlp_PT}(2,hlp,:),1,[])./reshape(E_pert_save_PT{hlp_PT}(2,hlp,:),1,[]),reshape(drabi_fit_save_PT{hlp_PT}(2,hlp,:),1,[])./reshape(E_pert_save_PT{hlp_PT}(2,hlp,:),1,[]),'o','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:),'MarkerFaceColor',cc(hlp,:));
            case {5}
                errorbar(alpha_sweep{hlp_PT}/pi*180,reshape(rabi_fit_save_PT{hlp_PT}(2,hlp,:),1,[])./reshape(E_pert_save_PT{hlp_PT}(2,hlp,:),1,[]),reshape(drabi_fit_save_PT{hlp_PT}(2,hlp,:),1,[])./reshape(E_pert_save_PT{hlp_PT}(2,hlp,:),1,[]),'s','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:),'MarkerFaceColor',cc(hlp,:));
            case {4,6}
                errorbar(alpha_sweep{hlp_PT}/pi*180,reshape(rabi_fit_save_PT{hlp_PT}(2,hlp,:),1,[])./reshape(E_pert_save_PT{hlp_PT}(2,hlp,:),1,[]),reshape(drabi_fit_save_PT{hlp_PT}(2,hlp,:),1,[])./reshape(E_pert_save_PT{hlp_PT}(2,hlp,:),1,[]),'s','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:));
            case {7,8}
                errorbar(alpha_sweep{hlp_PT}/pi*180,reshape(rabi_fit_save_PT{hlp_PT}(2,hlp,:),1,[])./reshape(E_pert_save_PT{hlp_PT}(2,hlp,:),1,[]),reshape(drabi_fit_save_PT{hlp_PT}(2,hlp,:),1,[])./reshape(E_pert_save_PT{hlp_PT}(2,hlp,:),1,[]),'^','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:));
            case {9}
                errorbar(alpha_sweep{hlp_PT}/pi*180,reshape(rabi_fit_save_PT{hlp_PT}(2,hlp,:),1,[])./reshape(E_pert_save_PT{hlp_PT}(2,hlp,:),1,[]),reshape(drabi_fit_save_PT{hlp_PT}(2,hlp,:),1,[])./reshape(E_pert_save_PT{hlp_PT}(2,hlp,:),1,[]),'^','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:),'MarkerFaceColor',cc(hlp,:));
            otherwise
                errorbar(alpha_sweep{hlp_PT}/pi*180,reshape(rabi_fit_save_PT{hlp_PT}(2,hlp,:),1,[])./reshape(E_pert_save_PT{hlp_PT}(2,hlp,:),1,[]),reshape(drabi_fit_save_PT{hlp_PT}(2,hlp,:),1,[])./reshape(E_pert_save_PT{hlp_PT}(2,hlp,:),1,[]),'o','MarkerSize',12,'LineWidth',1.5,'Color',cc(hlp,:));
        end
    end
    ymax=(max(max(rabi_estimate_save(2,:,:))))+0.25;
    ylim([0,ymax])
    xlim([0,90])
    xticks(0:15:90)
    xlabel('\alpha(\circ)','FontSize',18)
    title('DQ','FontSize',20)
    hold off
    set(gca,'FontSize',18)

    title('DQ','FontSize',20)
    hold off
    set(gca,'FontSize',18)
    
%     saveas(gcf,['DQ_',name_PT{hlp_PT}],'epsc')
    % plot g components
    
    figure
    cc_g=colormap(piratepal(9,1));
    hold on
    for hlp=1:6
        plot(alpha_fit/pi*180,g_munu_theory(hlp,:),'LineWidth',1.5,'Color',cc_g(hlp,:))
    end
    
    for hlp=1:6
        errorbar(alpha_sweep{hlp_PT}/pi*180,g_munu_PT{hlp_PT}(hlp,:),g_munu_err_PT{hlp_PT}(hlp,:),'o','MarkerSize',12,'LineWidth',1.5,'HandleVisibility','off','Color',cc_g(hlp,:))
    end
    hold off
    xlabel('\alpha(\circ)','FontSize',18)
    title('g_{\mu\nu}','FontSize',20)
    % legend('g_{\alpha\alpha}','g_{\phi+\phi+}','g_{\phi-\phi-}','g_{\alpha\phi+}','g_{\alpha\phi-}','g_{\phi+\phi-}','FontSize',12);
    set(gca,'FontSize',18)
    xlim([0,90])
    xticks(0:15:90)
    % ylim([0,6])
    
    lgd=legend('g_{\alpha\alpha}','g_{\beta\beta}','g_{\phi\phi}','g_{\alpha\beta}','g_{\alpha\phi}','g_{\beta\phi}');
    lgd.FontSize=18;
    
%     saveas(gcf,['metric_tensor_',name_PT{hlp_PT}],'epsc')
    %% calculate Berry curvature and the observable
    H_berry_curvature=zeros(1,length(alpha_sweep{hlp_PT}));
    H_berry_err=H_berry_curvature;
    
    for hlp_alpha=1:length(alpha_sweep{hlp_PT})
        [H_berry_curvature(hlp_alpha),H_berry_err(hlp_alpha)]=Herror_berry_curvature(g_munu_PT{hlp_PT}(:,hlp_alpha),g_munu_err_PT{hlp_PT}(:,hlp_alpha));
    end
    
    % Topological charge
    [QT,QT_err] = Herror_charge(alpha_sweep{hlp_PT},H_berry_curvature, H_berry_err);
    
    disp(['Q_T=',num2str(QT),' +- ',num2str(QT_err)])
end


function [rabi_estimate,g_munu_save,H_theory] = QGT_rabi_freq_detune_theory(alpha_sweep,beta,phi,rd,r,flag)

rabi_estimate = zeros(2,length(flag),length(alpha_sweep));
dpert=0.005;
g_out = zeros(length(flag),length(alpha_sweep));
g_munu_save = zeros(6,length(alpha_sweep));
H_theory = zeros(1,length(alpha_sweep));

for hlp=1:length(alpha_sweep)
    alpha=alpha_sweep(hlp);
    H0=Ham_detune(alpha,beta,phi,rd);
    [V,E]=eig(H0);
    [E,idx]=sort(diag(E));
    
    V=V(:,idx);
    um=V(:,1)/sqrt(V(:,1)'*V(:,1));
    u0=V(:,2)/sqrt(V(:,2)'*V(:,2));
    up=V(:,3)/sqrt(V(:,3)'*V(:,3));
    
    for hlp_flg=1:length(flag)
        flg_tmp=flag{hlp_flg};
        
        alpha_oscillate=flg_tmp(1);
        beta_oscillate=flg_tmp(2);
        phi_oscillate=flg_tmp(3);
        
        for hlp_omega=1:2 % omega
            dH=Ham_detune(alpha+dpert*alpha_oscillate,beta+dpert*beta_oscillate,phi+dpert*phi_oscillate,rd);
            dH=(dH-H0)/dpert;
            
            switch hlp_omega
                case 1
                    rabi_estimate(hlp_omega,hlp_flg,hlp)=abs(um'*dH*u0*r);
                    g_out(hlp_flg,hlp)=g_out(hlp_flg,hlp) + ...
                        abs(um'*dH*u0*r)^2/((E(1)-E(2))*r)^2;
                case 2
                    rabi_estimate(hlp_omega,hlp_flg,hlp)=abs(um'*dH*up*r);
                    g_out(hlp_flg,hlp)=g_out(hlp_flg,hlp) + ...
                        (rabi_estimate(hlp_omega,hlp_flg,hlp))^2/((E(1)-E(3))*r)^2;
            end
        end
    end
    
    g_munu_save(1,hlp)=g_out(1,hlp);
    g_munu_save(2,hlp)=g_out(2,hlp);
    g_munu_save(3,hlp)=g_out(3,hlp);
    g_munu_save(4,hlp)=(g_out(4,hlp)-g_out(5,hlp))/4;
    g_munu_save(5,hlp)=(g_out(6,hlp)-g_out(7,hlp))/4;
    g_munu_save(6,hlp)=(g_out(8,hlp)-g_out(9,hlp))/4;
    
    g_munu_tmp=[0,g_munu_save(4,hlp),g_munu_save(5,hlp);0,0,g_munu_save(6,hlp);0,0,0];
    g_munu_tmp=g_munu_tmp+g_munu_tmp.'+diag([g_munu_save(1,hlp),g_munu_save(2,hlp),g_munu_save(3,hlp)]);
    H_theory(hlp)=real(4*sqrt(det(g_munu_tmp)));
    
end
end

function Ham_out = Ham_detune(alpha,beta,phi,rd) %
Ham_out=[0,cos(alpha)*exp(-1i*beta),0;cos(alpha)*exp(1i*beta),0,sin(alpha)*exp(1i*phi);...
    0,sin(alpha)*exp(-1i*phi),0]+diag([rd,0,-rd]);
end



